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ABSTRACT 


This  thesis  presents  a  fully  coupled,  quasi-3D  analysis  of  rotor  blade  flutter  that 
can  accommodate  forward  flight  conditions.  The  rotor  blade  is  modeled  as  a  uniform 
beam,  taking  the  average  characteristics  of  a  real  blade  between  20%  and  90%  of  its 
length.  Applying  Rayleigh’s  method,  the  first  few  bending  and  torsion  normal  mode 
shapes  and  natural  frequencies  are  detennined,  and  then  adjusted  for  the  rotating  case. 
With  this  data,  force  and  moment  equations  of  motion  are  developed  using  Lagrange’s 
equation  along  with  a  normal  mode  analysis.  Theodorsen  coefficients  are  calculated  over 
a  range  of  forward  velocities  (input  as  reduced  frequencies)  for  a  specified  number  of 
elements  along  the  blade  model.  Incorporating  these  coefficients  into  the  equations  of 
motion,  a  square  matrix  is  generated  from  which  complex  eigenvalues  can  be  derived. 
These  eigenvalues  provide  the  aeroelastic  natural  frequencies  and  damping  coefficients 
for  each  coupled  mode.  The  forward  velocity  at  which  one  of  the  modes  produces  a 
positive  damping  coefficient  gives  the  value  of  reduced  frequency  for  the  flutter  point. 
The  resulting  forward  speed  and  blade  tip  speed  can  then  be  determined. 
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I.  INTRODUCTION 


A.  INTRODUCTION 

The  conventional  method  for  designing  a  rotor  blade  to  be  free  of  flutter 
throughout  the  helicopter’s  flight  regime  is  to  design  the  blade  so  that  aerodynamic  center 
(AC),  elastic  axis  (EA)  and  center  of  gravity  (CG)  are  coincident  and  located  at  the 
quarter-chord.  This  practice  pays  off  by  decoupling  the  equations  used  in  two- 
dimensional  unsteady  aerodynamic  theory.  While  this  assures  freedom  from  flutter,  it 
adds  constraints  on  rotor  blade  design  not  usually  followed  in  fixed  wing  design.  A  blade 
designed  with  the  CG,  EA  and  AC  coincident  at  the  quarter  chord  is  heavier  than  one  free 
of  that  restriction.  It  also,  if  strictly  followed,  rules  out  use  of  a  flap  that  causes  the  AC  to 
move  with  flap  angle.  In  addition,  it  restricts  use  of  camber  that  moves  the  AC  aft  [Ref. 
1]. 

Loewy’s  [Ref.  2]  2-D  unsteady  aerodynamic  theory,  as  amended  by  Jones  and 
Rao  [Ref.  3]  and  Hammond  [Ref.  4],  provides  a  useful  tool  for  examining  blade  flutter  in 
a  hover.  Extension  of  their  work  to  a  helicopter  in  forward  flight  presents  a  fonnidable 
mathematical  challenge,  and  thus,  at  present,  there  is  no  accepted  theory  to  completely 
analyze  blade  flutter  in  forward  flight.  Currently,  to  meet  forward  flight  blade  flutter 
requirements  the  rotorcraft  manufacturer  must  rely  on:  (1)  a  quasi-fixed  wing  blade 
flutter  analysis,  which  does  not  account  for  the  unsteady  contribution  of  the  wake  below 
the  rotor;  and  (2)  costly  rotor  whirl  tests  at  normal  and  over-speed  conditions  that,  while 
providing  information  in  regard  to  blade  flutter,  do  not  accurately  simulate  either  blade 
dynamics  or  unsteady  aerodynamics  in  forward  flight. 

However,  closer  examination  of  the  problem  reveals  that  it  is  possible  to  make 
several  simplifying  assumptions  that  make  the  forward  flight  flutter  problem  tractable.  In 
particular,  it  is  assumed  that  at  the  onset  of  flutter  oscillations  begin  to  build  up  prior  to 
the  blade  reaching  a  critical  azimuth  position,  then  decay  as  the  blade  moves  beyond  this 
point.  Shipman  and  Wood  [Ref.  5]  provide  the  two-dimensional  basis  for  this  three- 
dimensional  analysis.  This  thesis  will  follow  the  analysis  of  Shipman  and  Wood,  using 
Theodorsen’s  lift  deficiency  function  (2-D  approach),  and  further  expand  the  theory  to 
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encompass  full  blade-length  dynamics  (3-D  effects).  A  numerical  example  using  UH-60 
blade  data  is  included  as  a  means  of  illustration  and  validation. 

B.  BACKGROUND 

“ Main  and  tail  rotor  blades  of  the  A.T.H.  have  been 
designed  so  that  center  of  gravity,  elastic  axis,  and  aerodynamic 
center  are  coincident.  Also,  the  control  system  for  the  main  rotor 
is  stiff  with  high  internal  damping.  No  main  or  tail  rotor  blade 
flutter  has  been  experienced  with  earlier  model  helicopters 
possessing  these  design  features. 

Main  and  tail  rotor  blades  for  the  HSS-2,  which  are  the 
same  as  those  of  the  A.T.H. ,  have  been  installed  on  Sikorsky  whirl 
stands,  and  tested  at  maximum  design-limit  speeds.  Main  rotor 
blades  were  tested  for  power-on  and  power-off  conditions.  Tail 
rotor  test  conditions  were  power-on  and  power-off.  Observation  of 
blades  during  these  tests  indicated  no  flutter  or  divergence  at 
maximum  operating  conditions.  ”  [Ref.  6] 

This  statement,  from  the  1960  Sikorsky  Report  No.  50131  for  the  Advanced 
Tactical  Helicopter,  is  representative  of  the  current  methods  for  ascertaining  freedom 
from  flutter  for  new  helicopter  acquisitions.  Theoretical  calculations  consist  of  a  2-D 
flutter  analysis  of  the  blade  tip  as  if  it  was  a  fixed  wing,  and  whirl  stand  tests  involve 
significant  expense.  Much  effort  has  been  expended  in  the  last  century  to  improve  flutter 
prediction  and  minimize  the  cost  in  time  and  money  spent  on  whirl  stand  tests.  The 
following  paragraphs  are  a  review  of  some  of  the  significant  work  that  has  advanced  the 
study  of  flutter  and  provided  methods  of  solving  the  flutter  problem. 

1,  Theodorsen  and  Loewy 

In  1936,  Theodorsen  [Ref.  7]  made  the  flutter  problem  somewhat  tractable  with 
certain  simplifying  assumptions.  He  considered  a  wing  of  infinite  aspect  ratio, 
encountering  small  oscillations  at  a  constant  velocity  through  an  incompressible,  non- 
viscous  fluid.  By  using  these  assumptions,  the  forces  and  moments  could  be  detennined 
from  2-D  potential  flow  theory.  In  addition,  instead  of  using  the  actual  mass  and 
geometry  distributions,  he  assumed  that  the  results  could  be  conservatively  obtained  for  a 
unit  span  of  the  wing.  Using  the  nomenclature  found  in  Scanlan  and  Rosenbaum  [Ref. 
8],  and  excluding  tenns  associated  with  an  aileron,  the  unsteady  force  and  moment  per 
unit  span  are  given  by 
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The  term  C(k)  is  Theodorsen’s  well-known  lift  deficiency  function  defined  by 
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where  H|2)  =  Jn  -  iYn  is  the  Hankel  function  of  the  second  kind  of  order  n  evaluated  at 

reduced  frequency  k  =  cob/V.  The  force  and  moment  equations  take  into  account  both 
pitch  (torsion)  and  plunge  (flapping)  motion,  and  include  circulatory  and  non-circulatory 
terms.  It  is  in  the  circulatory  terms  that  the  lift  deficiency  function  is  defined. 

Inherent  in  Theodorsen’s  2-D  thin  airfoil  theory  is  the  assumption  that  the  wake  is 

carried  downstream  to  infinity  by  the  free-stream  airflow  [Ref.  9].  For  work  in  rotor 
blade  flutter,  this  assumption  presents  problems,  since  the  rotor  blade  sections  may 
encounter  the  shed  wake  from  previous  blades  as  well  as  the  reference  airfoil.  Loewy 
[Ref.  2]  recognized  this  problem  and  developed  a  2-D  model  for  the  hover  case  that 
accounted  for  the  effects  of  previously  shed  wakes.  The  wakes  were  modeled  as  a  series 
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of  planar  2-D  vortex  sheets  separated  by  a  distance  ‘h’  that  was  dependent  on  the  induced 
velocity  through  the  rotor  disk  and  the  number  of  rotor  blades.  Loewy  then  developed  a 
modified  lift  deficiency  function  that  could  be  used  in  conjunction  with  the  force  and 
moment  equations.  This  modified  lift  deficiency  function  is  defined  as 

C'(k  m  h)  = _ H\2\k)  +  2J\ (kW_ (k,m, h) _ 

where 

W(k,m,h )  =  —h — - - 

eTe^m)_1 

and  is  evaluated  at  reduced  frequency  (k),  wake  spacing  (h  =  InV/bQQ.)  and  frequency 
ratio  (m  =  CO /£2).  Note  that  as  h  — >  °°  ,  W  — >  0  and  C'(k)  — >  C(k) . 

2.  Shipman  and  Wood 

In  1971,  Shipman  and  Wood  [Ref.  5]  examined  the  unsteady  aerodynamics  acting 
on  an  advancing  blade  in  steady-state  level  flight.  The  basic  assumptions  of  their  work 
are  as  follows: 

1.  Two-dimensional,  inviscid,  incompressible  potential  flow. 

2.  Respective  layers  of  the  wake  are  two-dimensionalized  and  treated  as  parallel 
horizontal  sheets. 

3.  In  forward  flight,  each  blade  of  the  rotor  will  respond  in  the  same  manner  as 
every  other  blade. 

4.  The  most  critical  azimuth  position  of  the  blade  in  forward  flight  for  the  onset 
of  flutter  is  at  y/=  90  . 

5.  At  the  onset  of  blade  flutter,  oscillations  will  begin  to  build  up  prior  to  the 
blade  reaching  the  critical  azimuth  position,  and  these  oscillations  will  decay 
as  the  blade  moves  beyond  the  critical  azimuth  position. 

At  a  specified  radial  location  r  on  the  blade,  the  local  tangential  velocity  is  given 


Ut  (r)  =  £2r  +  V  sin\|/ 

If  it  is  assumed  that  the  flutter  speed  for  this  blade  segment  is  such  that 

UFL(r)<ftr  +  V 
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then  during  each  blade  revolution  the  blade  segment  at  radius  ‘r’  will  experience 
velocities  which  will  increase  to  the  flutter  speed  and  beyond,  then  return  through  the 
flutter  point  to  lower  airspeeds.  Extending  this  concept  to  both  blade  azimuth  position 
and  radial  position,  one  can  observe  that  the  blade  tangential  velocity  at  a  given  radial 
position  will  exceed  the  flutter  speed  in  some  region  of  rotor  azimuth  position  if: 

V sin\|/  >  UFL  (r)  -  fir 
This  flutter  region  is  illustrated  in  Fig.  1 . 


Figure  1.  Unstable  Region  on  Advancing  Blade 

It  should  be  noted  that  all  points  within  the  shaded  region  of  Fig.  1  experience 
negative  damping.  This  negative  damping  when  combined  with  transient  oscillations  in 
the  blade  and  will  cause  blade  motion  to  grow.  Also,  in  the  region  \) /  <  Till  -  A\|/FL,  the 
damping  will  decrease  as  \|/  approaches  (7t/2  -  A\|/FL),  whereas  in  the  region  \\f  >  Till  + 
A\|/FL,  damping  will  be  positive  and  will  increase  such  that  blade  transients  tend  to  die 
out. 

Consider  the  effect  of  this  variation  in  damping  on  an  outboard  portion  of  the 
advancing  blade.  It  is  expected  that  damping  would  decrease  as  the  blade  approaches  \| / 

o 

=  90  and  the  amplitude  of  oscillations  would  build  up.  Conversely,  as  the  blade 

o 

advanced  beyond  \|/  =  90  ,  damping  would  increase,  and  there  would  be  a  corresponding 
decrease  in  blade  vibratory  amplitude.  This  build-up  and  decay  of  blade  amplitude  would 
result  in  a  distribution  of  shed  vorticity  as  shown  by  Fig.  2.  Here,  we  observe  that  time- 
wise  variations  in  amplitude  of  blade  vibrations  have  resulted  in  space- wise  variations  in 
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shed  vorticity.  Since  we  have  assumed  steady-state  flight,  each  blade  would  shed  similar 
segments  of  vorticity  for  each  revolution.  These  vortex  segments  constitute  the  wake  that 
will  be  treated  in  this  analysis. 


Figure  2.  Distribution  of  Shed  Vorticity  in  Unstable  Region 

Based  on  the  foregoing,  the  bound  vorticity  on  the  airfoil  can  be  expressed  as  the 
product  of  a  function  of  chord-wise  position,  a  decay  function,  and  a  harmonic  function 
of  time1.  We  write  the  incremental  bound  vorticity  as 

Ya=Ya(^feo)e'>^) 

where  f(^0)  is  an  assumed  decay  function  centered  about  =  0.  The  limiting  case  of 
constant-strength  shed  vorticity  such  as  considered  by  Theodorsen  [Ref.  7]  and  Loewy 
[Ref.  1]  for  their  analyses,  is  simply  achieved  by  taking  f(q0 )  =  1. 

When  the  inflow  velocity  through  the  rotor  is  small,  the  shed  vorticity  remains 
close  to  the  rotor  and  the  wakes  shed  from  each  blade  during  several  previous  passes  as 

well  as  the  present  pass  must  be  considered.  The  build-up  and  decay  of  vorticity  occurs 

° 

within  a  double  azimuth  angle  on  either  side  of  \|/  =  90  .  The  solid  lines  of  Fig.  3 
indicate  this  region  of  the  wake.  In  this  region  the  azimuth  angle  between  a  shed  vortex 
filament  and  the  reference  blade  may  be  ignored.  The  tip  does  not  move  very  far  from 
the  vertical  plane  shown  in  Fig.  3  and  so  its  path  may  be  taken  to  lie  in  this  plane. 


1  W.  P.  Jones  [Ref.  10]  first  treated  the  case  of  an  oscillating  airfoil  where  the  strength  of  the  airfoil’s 
position  was  allowed  to  grow  or  decay  exponentially  with  time.  As  the  buildup  or  decay  rate  approached 
zero,  Jones’  lift  deficiency  function  approached  that  of  Theodorsen  [Ref.  6], 
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Figure  3.  Development  of  Skewed  Helical  Wake 

Combining  the  vorticity  segments  given  in  Fig.  2,  the  resulting  wake  pattern  is 
shown  in  Fig.  4.  With  the  mathematical  model  defined,  the  problem  now  is  to  determine 
the  pressure  difference  across  the  airfoil  due  to  the  vorticity  shed  in  the  wake,  and 
consequently  to  determine  the  unsteady  lift  and  moment  acting  on  the  airfoil. 


Figure  4.  2-D  Wake  Model  for  Forward  Flight 
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Shipman  and  Wood  developed  their  equation  for  the  unsteady  forces  in  forward 
flight  that  are  analogous  to  Theodorsen  and  Loewy,  but  modified  the  lift  deficiency 
function  to  account  for  the  helicopter’s  forward  speed  (advance  ratio)  and  the  build-up 
and  decay  function  associated  with  the  advancing  blade  illustrated  in  Fig.  2.  The  forward 
flight  lift  deficiency  function  is  defined  by 

_  ,,  7  .  H\2\k)  +  AFAk)  +  AFAk)  +  2JAk)  \W(k,h,m)  + AW(k,h,m)] 

f  (  /'  Li  711  |  —  _ _ _ _ _ 

1  H^(k)  +  iHf\k)  +  AFfk)  +  2{Jl(k)  +  iJ0(k))  [W(k,h,m)  +  AW(k,h,m)\ 


where 


AF2(k) 

AFAk) 

AFAk) 


2 

n 

2 

k 

2 

7T 


,  1  df  ye-'ky 

ik  dy  2  _1 


dy 


\  .  1  df 

ik  dy 


y  + 1 
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-iky 


df  1  d  f 


dy  ik  dy2 


y-\ 


dy 

iky  dy 


AW(k,h,m) 
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—knh  iknm 

e  e 


n= 1 


X 


/  (—nm  -  inh )  —  1  — 


1  df  (—nm  -  inh  ) 
ik  dy 


where  f(y)  is  the  assumed  decay  function.  Note  that  as  h  — ■>  °°  ,  and  5  — »  °°  , 
Cfk)  C(k) ,  and  when  V  -»  0,  Cfk)  -»  C'(fc) . 
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II.  THEORETICAL  DEVELOPMENT 


A.  A  FLUTTER  THEORY  FOR  FORWARD  FLIGHT 


Shipman  and  Wood,  using  the  theory  described  in  the  previous  section  derived  an 
equation  analogous  to  Theodorsen  and  Loewy  that  accounted  for  the  effects  of  forward 
flight  and  a  shed  skewed  helical  wake.  However,  evaluation  of  some  of  the  integrals 
required  for  the  AF  and  AW  terms  is  extremely  difficult.  If  the  view  is  taken  that  the 
Theodorsen  method,  which  neglects  the  effect  of  the  shed  wakes,  is  just  a  first-order 
approximation  to  the  rotary-wing  flutter  problem,  then  use  of  the  Theodorsen  lift 
deficiency  function  should  yield  suitable  first-order  results.  Therefore,  this  thesis  will 
rely  on  Theodorsen’ s  lift  deficiency  function,  while  following  the  methods  of  Shipman 
and  Wood  in  modeling  forward  flight. 

1.  Mode  Shapes  and  Natural  Frequencies 

Blade  bending  frequencies  and  mode  shapes  used  in  this  analysis  were  determined 
by  a  simplified  method  by  assuming  a  unifonn  stiffness  and  mass  distribution  along  the 
blade.  A  more  exact  and  detailed  analysis  would  be  required  to  account  for  such  details 
such  as  local  changes  in  stiffness  and  mass  distribution  due  to  blade  features  such  as 
doublers  near  the  blade  root  and  outboard  blade  balance  weights  near  the  tip.  While  these 
details  would  be  desired  for  an  actual  blade  design,  they  can  be  viewed  as  second  order 
effects  and  not  necessary  for  the  first  order  flutter  analysis. 

For  the  simplified  model,  the  geometric  and  inertial  properties  of  the  subject  blade 
are  averaged  between  20%  and  90%  of  its  length.  The  Fourier-based  solution  of  the 
pinned-free  uniform  beam  from  Young  and  Felgar  [Ref.  1 1]  are  applied  to  determine  the 
non-rotating  mode  shapes: 


yn  =  ^-(cosh(r)  +  cos(r)  -  An  (sinh(r)  +  sin(r))) 


where  An  =  1.000777,  1.000001  and  1.0  for  the  first  three  bending  modes.  Then  the 
method  given  by  Den  Hartog  [Ref.  12]  for  determining  non-rotating  natural  frequencies 
is  applied: 
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CO  =  a 


where  an  =  15.4182,  49.9649  and  104.2477  for  the  first  three  bending  inodes,  and  |U  is  the 
mass  per  unit  length. 

Given  the  natural  frequencies,  the  rotational  velocity  of  the  rotor  head,  and  the 
non-rotating  mode  shapes,  Yntema  [Ref.  13]  is  called  on  to  supply  the  rotating  natural 
frequencies.  Figures  5  through  7,  taken  from  Yntema’s  report,  compares  rotating  and 
non-rotating  mode  shapes  for  the  first  three  bending  modes  of  a  pinned-free  beam  to 
validate  the  assumption  that  the  non-rotating  mode  shape  is  a  close  approximation  of  the 
rotating  mode  shape. 
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1.0 
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Yntema  notes  that  an  exact  value  for  the  nth  bending  frequency  of  a  beam  rotating 
at  any  rotational  speed,  Q,  can  be  found  if  the  nth  natural  bending  mode  shape  is  known 
for  this  value  of  rotational  speed.  He  obtains  his  frequency  equation  by  equating  the 
kinetic  energy  at  zero  displacement  to  the  potential  energy  of  both  the  bending  and 
centrifugal  forces  at  maximum  displacement  for  vibration  perpendicular  to  the  plane  of 
rotation. 


8,J 

oEI>V  dx  ,  J 

•1  ,ry 

ij.y,  * 

f'my;  dx  i 

[  mYn  dx 

where  n  refers  to  mode  of  vibration  and 


(1) 


Ti=j(rl  +  e)m  drj 


He  goes  on  to  point  out  that  while  the  rotating  mode  shape  is  unknown,  a  close 
approximation  to  the  rotating  natural  frequency  can  be  obtained  by  making  use  of 
Rayleigh ’s  Principle,  and  using  the  non-rotating  mode  shape  in  equation  (1).  The  report 
states  that  the  non-rotating  mode  shape  is  consistent  with  the  constraints  of  the  system  (in 
this  case  a  pinned-free  beam).  If  the  nth  mode  of  the  non-rotating  mode  shape,  yn,  is 
substituted  into  equation  (1),  the  first  term  becomes  exactly  the  square  of  the  bending 
frequency  of  the  non-rotating  beam.  By  denoting  the  ratio  in  the  second  term  by  Kn,  a 
Southwell  coefficient,  the  form  of  the  frequency  equation  becomes: 

(D2  =(D2  +KQ2  (2) 

R  NR 

n  n 

To  account  for  blade  offset,  e,  subdivide  Kninto  two  independent  parts: 

Kn  =  K0  +  Kj  e  (3) 

11  ln 


where  K0n  is  referred  to  as  the  zero-offset  Southwell  coefficient  and  Kin  is  referred  to  as 
the  offset-correction  factor  for  the  Southwell  coefficient.  As  is  frequently  done,  it  is 
convenient  to  write  the  square  of  the  non-rotating  frequency  in  terms  of  a  non-rotating- 
beam  frequency  coefficient,  an,  and  the  mass  and  stiffness  of  the  beam  as: 


EL 


m0l 


(4) 
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Combining  equations  (2),  (3)  and  (4)  yield: 


to; 


K  ~  ®nr„  +  (ko„  +Klne)^2 


Yntema’s  report  [Ref.  13]  gives  charts  that  provide  an,  K0n  and  Ki„  which,  in 
conjunction  with  the  mass  and  stiffness  of  the  beam  at  the  root,  the  length  of  the  beam, 
the  hinge  offset,  and  the  rotational  speed,  permit  rapid  estimation  of  the  first  three 
bending  frequencies  of  rotating  beams  with  hinged  or  cantilevered  root-end  support.  As 
previously  noted,  once  the  frequencies  have  been  found,  the  rotating  beam  mode  shapes 
can  then  be  approximated  by  non-rotating  mode  shapes  which  are  defined  by  the  Fourier- 
based  solutions  contained  in  Young  and  Felgar  [Ref.  11]. 

The  next  problem  is  to  determine  the  mode  shapes  and  frequencies  for  the 
uniform  beam  blade  model  in  the  torsional  case.  For  this,  Den  Hartog  [Ref.  13]  provides 
a  relatively  simple  solution.  Given  that  the  torsional  stiffness  along  the  beam  is  constant, 
the  torsional  mode  shapes  are  given  by 


y„  =  sm 


n  - 


1 


7tr 


R 


The  non-rotating  torsional  natural  frequencies  are  given  by 


co, 


NRn 


(  n 

" 

n  — 

7C  - 

l  2; 

i 

GJ 

fji2 


Bramwell  [Ref.  14],  gives  the  exact  solution  for  rotating  natural  frequencies  in  the 
torsional  case  as 


®R„  ®NR„  ^ 


2.  Normal  Mode  Analysis 

With  the  rotor  blade  free  vibration  problem  solved,  the  bending  deflection  may  be 
defined  as 
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h(X’t)  =  Zfn(X)Cln(t) 

n-1 

where  fn(x)  is  the  characteristic  function  (mode  shape)  for  the  nth  vertical  bending  mode 
of  the  rotor  blade.  The  quantities  qn(t)  can  be  considered  as  weighting  functions  for  each 
mode  that  contributes  to  the  deflection.  They  are  called  the  normal  coordinates  since 
they  can  be  shown  to  reduce  the  kinetic  and  potential  energy  expressions  to  sums  of 
squares  of  the  coordinates  with  no  cross  product  tenns. 

The  corresponding  torsional  deflection  of  the  rotor  blade  can  be  written  in  terms 
of  the  blade  torsion  modes  as: 

a(x’t)  =  ZFn(x)cln(t) 

n=l 

where  Fn(x)  is  the  characteristic  function  of  the  nth  torsional  mode  of  the  rotor  blade  and 
qn(t)  is  the  corresponding  nonnal  coordinate.  Consider  the  five-degree  of  freedom  (DOF) 
case  with  three  bending  modes  and  two  torsion  modes.  The  bending  and  torsional 
deflections  can  be  written  as: 

h(x,t)  =  hj  (t)fj  (x)  +  h2  (t)f2  (x)  +  h3  (t)f3  (x) 
and 

a(x,t)  =  ai(t)Fi(x)  +  a2(t)F2(x) 

where 

fi(x)  =  1st  vertical  bending  mode 
f>(x)  =  2nd  vertical  bending  mode 
f?(x)  =  3ld  vertical  bending  mode 
Fi(x)  =  1st  torsion  mode 
F2(x)  =  2nd  torsion  mode 

Now  all  the  tools  are  in  place  to  begin  the  detailed  flutter  analysis,  which  is  the 
combination  of  Lagrange’s  equation  with  Theodorsen’s  work  to  provide  the  equations  of 
motion. 
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3.  Lagrange  and  The  Equations  of  Motion 

Lagrange’s  equation  is  given  as: 


df  3T  ^ 
dt^qny 


3T  |  3U  |  3D 

dqn  dqn  dqn 


where  T  =  kinetic  energy,  U  =  potential  energy,  D  =  dissipation  function,  and  Qn  = 
generalized  force.  For  the  five  DOF  case, 

T  =  t M,hf  +  iM,h=  +  +  tl„a?  +  a; 

+S„ih,a, +S^h,a2  +  S„!ih2a, 

+^“22^2^2  +  S0,3ih3(X1  +Sa^h3(X2 


U  =  -Mjtog  hf  +  -M2(diy2  +  -M3<h32 

+  |1«1<«f  +|Ia2<«2 


D  =  +  M,gh2<h^  +  M3ghi  <  h' 


2co 


2co 


2co 


|  Iq,g „,<«!  [  Iq2g«2<«2 


2co 


2co 


and  the  generalized  forces  including  aerodynamic  terms  are  defined  as: 


Qh,  —  tip co  ( Ajjh,  +  A12h2  +  A13h3  +  Ai40Cj  +  A15cx2 ) 
Qh,  —  ttpco  (Ajjhj  +  A22h2  +  A23h3  +  A240Cj  +  A25cx2) 

Qh3  —  tipco  (A31h,  +  A32h2  +  A33h3  +  A340Cj  +  A35cx2) 

Qa,  —  ttpco  (A41h,  +  A42h2  +  A43h3  +  A440Cj  +  A45oc2) 

Q«2  —  ttpco  ( A51h,  +  A52h2  +  A53h3  +  A54oc,  +  A55oc2 ) 


The  expressions  for  aerodynamic  terms  that  couple  the  modes  together  and  incorporate 
Theodorsen’s  lift  deficiency  function  are  given  as 
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(5) 


An  =  Jo'b2  [fi  (x)]2  Lhdx 

Ai2=|01b2f1(x)f2(x)Lhdx 

A13  =  J0b2fl  (X)f3(X)Lhdx 


A 


A 


14 


15 


A 


A 


24 


25 


A 


34 


35 


j0b3fi(x)F1(x) 

La- 

'l  ) 

v2  J 

|0b3fi  (x)F2  (x) 

La" 

f  1 

u  J 

A21  =  (X)fl  (X)MX  =  A12 


A22  =Job2[f2(x)]"Lhdx 
A23=J01b2f2(x)f3(x)Lhdx 


j0b3f2(x)F1(x) 

K~ 

v2  J 

j0b3f2(x)F2(x) 

L." 

U  J 

A31  =  J0b2f3  (x)fi  (x)Lhdx  =  A13 
A32  =  (X)f2  (X)Lhdx  =  A23 


A33  =  Job2  [f3  (X)]2  Lhdx 


J0*  b3f3  (X)F,  (x) 

L»  “I 

fl  ) 

1 2  J 

J0b3f3  (x)F2  (x) 

La  “I 

f  1 

U  J 

dx 


dx 


dx 


dx 


dx 


dx 
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A4i  =£b3F1(x)f1(x) 

A42  =Jl'b3F1(x)f2(x) 


Mh- 


M, 


f  1  > 

— +  a 

v2  j 


L„ 


f\  ^ 

—  +  a 

v  2  , 


L„ 


dx 


dx 


M, 


A44=j0b4[Fi(x)] 

A45  =j01b4F1(x)F2(x) 


A43={ob3F1(x)f3(x) 

f\  ' 


f\  ^ 

— +  a 

v2  j 


L, 


dx 


M„  - 


—  +  a 

v  2  j 


(La  +Mh)  + 


1 

—  +  a 

v2  j 


U 


M„ 


f  1  ^ 

—  +  a 

V2  j 


(L«  +  Mh)  + 


1 

— +  a 

V  2  j 


L, 


dx 


dx 


Asi  =J01b3F2(x)f1(x) 

A52  =|01b3F2(x)f2(x) 


A53  =|0lb3F2(x)f3(x) 


M. 


Mh- 


Mh- 


— +  a 

v2  j 


L, 


\ 

v2  J 

( 1 

U  J 


L„ 


L, 


dx 


dx 


dx 


A54  =  J"q  b4F2  ( x )  Fj  ( x ) 


M„  - 


—  +  a 

v2  j 


(La  +Mh)  + 


f  1  ^ 

—  +  a 

V2  j 


L„ 


dx  =  A 


45 


A55=j0lb4[F2(x)]J 


M, 


f\  ^ 

—  +  a 

V2  j 


(La  +Mh)  + 


f\  > 

—  +  a 

V  2  j 


L. 


dx 


The  generalized  masses  of  the  three  bending  modes  and  two  torsion  modes  can  be 
written  as: 


Mr  =  Jom(x)[fi  (x)]2  dx 

M2  =  J*om(x)[f2  (x)]2  dx 

M3  =  Jom(X)[f3  (X)]  dX 
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Ia1  =j0Ia(X)[F1(X)]“dx 

la2  =j0'Ia(X)[F2(X)]“dx 

The  coupled  static  unbalance  terms  are  defined  as: 

sa„  =J01S„(x)f1(x)F1(x)dx 

Sa12  =j’Sa(x)f1(x)F2(x)dx 

Sa21  =J01S„(x)f2(x)F1(x)dx 

Sa22  =J0'Sa(x)f2(x)F2(x)dx 
Sa31  =j0'Sa(x)f3(x)F1(x)dx 
Sa32  =J01Sa(x)f3(x)F2(x)dx 

First  it  is  noted  that  the  kinetic  energy  equation  is  only  a  function  of  the  derivative 
of  the  generalized  displacement  ( hn  or  an ).  Thus,  Lagrange’s  equation  reduces  to: 

d  f  dT^  dU  |  3D  _ 
dt  ^  3qn  J  dqn  dqn 

Applying  Lagrange’s  equation  to  each  of  the  five  DOFs  yields  the  following  live 
equations: 


M|h,  +sa  a,  +s  a,  +  h,  + - -Qh 

CO 


M2h2  +S„iia,  +S„„«2  +  M2<h2  +  h2  -  Qh 

CO 


Mjhj  +  S„  ft,  +S  «2  +  M,<h,  +  M,<0|'-gt-  h,  =  Qh 
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+Sanhl  +Sa21h2  +Sa„h3  +Ial0)a,ai  + 


2  ,  ^a,  ®a,  §a. 


CO 


a.  =Qo 


!«,  +  Sa„  \  +  Sa„  h2  +  S  h3  + 1  <  a2  + 


^  u  .  ®a2  Sa2 

CO 


«2=Qa, 


If  simple  harmonic  motion  is  assumed,  that  is:  hn  =  -co2hn,  hn  =  icohn,  an  =  -co2an,  and 
an  =icoan,  and  the  expressions  for  Qhn  and  Q(/|i  are  substituted  into  the  equations  of 
motion,  the  results  are: 


rcpAn  +  M,  —  M,  (l  +  igh 


CO,, 


CO 


h1+(7tpA12)h2+(7tpA13)h3 


+ ( TUP  A14  +  San )  a,  +  (tip  A1S  +  Sai2 )  a2  =  0 


(7tpA21)h1  + 


7tpA22  +  M2-M2(l  +  igh 


CO,, 


CO 


h2  +  (ttpA23)h3 


+ 


(7tpA24  +  S„2]  )a,  +  (7tpA25  +  Sa22  )a2  =  0 


(7tpA31)h1  +(^pA32)h2  + 


7tpA33  +  M3-M3(l  +  ighi 


CO,, 


A2 


CO 


+ 


(7tpA34  +  S„31  )a,  +  (7tpA35  +  S„32  )a2  =  0 


(ttpA41  +  San  )h,  +  (7tpA42  +  Sa2i  )h2  +  (kPA43  +  Sa3]  )h3 


+ 


TtpA44  +  Iai  -Iai(l  +  ig0 


co„ 


CO 


a,  +  ;tpA45a2  =  0 


(7UpA51  +  Sa|2  )h,  +  (7tpA52  +  S„22  )h2  +  (7tpA53  +  S„32  )h3 


+;tpA54a1  + 


^PA55+Ia:  (l  +  ig0 


co„ 


CO 


a2  =0 
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The  five  equations  to  the  flutter  problem  can  be  written  in  matrix  fonn  as 


(6) 


[npAjj+M, 

-Mi(l+igh) 

7tpA2i 

^pA31 


*PA4i  +saii 


^PA5I  +\2 


npAn 

[rcpA^+M, 
-M2(l  +  ighJ 

7cpA32 

’«PA«+Stta 

^PA52+Stt22 


CO 


7tpAi3 

TtpA23 

[;ipA33+M3 
-M3(l+igJ 

’VAJJ+S«ta 


°V 

CO 


TtpAi4+Sa 
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[^44+^ 
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a, 

2 
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rcpAlS+Sa 
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KPA15+\2 


7TpA4, 
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[rcpAcs+Io, 
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4.  Eigenvalues 

Equation  (6)  is  a  complex  eigenvalue  problem.  In  order  to  solve  the  flutter 
problem,  it  is  convenient  to  rewrite  the  problem  in  the  form  ( A  -  IZ)X  =  0  by  arbitrarily 
letting 


( co„ 


V  ©  J 


(l  +  ig) 


which  results  in 


>\ 

l 

N 

CN 

A13 

A14 

a15 

^21 

A  22  —  Z 

A23 

A  24 

A25 

A31 

A32 

A33  -Z 

A34 

A35 

>1 

A  42 

A43 

A44  —  z 

A45 

A51 

A  52 

A53 

A  54 

A55  -Z 
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It  should  be  noted  that  the  coefficients  of  the  characteristic  equation  of  the 
(A-IZ)  matrix  (a  quintic  in  Z)  are  complex,  due  to  Theodorsen’s  lift  deficiency 

function,  and  thus  the  eigenvalues  will  be  complex.  The  coupled  frequency  of  oscillation 
(Cfl)  for  each  eigenvalue  can  be  found  from  the  real  part  of  Z,  since  the  first  torsional 
natural  frequency  is  already  known: 


CO;  =  ■ 

VRe(z) 

The  damping  coefficient  required  for  flutter  to  exist  (g)  for  each  eigenvalue  can 
be  found  from  the  imaginary  part  of  Z: 

r  \2 

g,  =  hn(Z) 

If  g  is  negative  for  the  reduced  frequency  chosen,  then  damping  must  be 
decreased  in  order  to  obtain  flutter.  Negative  values  of  g  represent  the  stable,  or  non¬ 
flutter,  condition.  If  g  is  positive,  then  damping  must  be  increased  to  be  at  the  flutter 
point.  Positive  values  of  g  represent  the  unstable,  or  flutter  condition.  When  a  plot  of  g 
is  made  against  1/k  (k  being  reduce  frequency),  there  will  be  five  curves  corresponding  to 
the  variation  of  each  eigenvalue  as  the  reduced  frequency  varies.  Some  of  these  curves 
will  have  only  values  of  g  that  are  negative.  These  are  the  non-critical  curves  that 
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co„ 


represent  stable  coupled  modes  and  do  not  influence  the  flutter  solution.  However,  at 
least  one  curve  will  start  with  a  negative  value  of  g  and  then  at  some  point  cross  the 
abscissa  (1/k)  to  a  positive  value  of  g.  This  curve  is  called  the  critical  curve,  and  the 
value  of  1/k  where  this  curve  crosses  the  abscissa  represents  the  critical  flutter  speed,  or 
flutter  point.  The  critical  flutter  speed  is  found  from  the  relationship: 


where  co  is  found  from  the  real  part  of  the  eigenvalue  relationship  described  above  for  the 
critical  curve  evaluated  at  the  reduced  frequency  that  corresponds  to  the  crossover  point 
(kcrit).  Once  the  unstable  mode  is  identified,  results  are  commonly  plotted  as  g  vs.  Ufl- 
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III.  RESULTS 


A.  UH-60  BLADE  EXAMPLE 

Using  the  data  accumulated  by  NASA  for  Sikorsky’s  UH-60,  a  sample  numerical 
problem  will  now  be  presented.  The  UH-60’s  blade  is  modeled  as  a  uniform  beam, 
incorporating  the  average  geometric  and  inertial  characteristics  of  the  blade  between  20% 
and  90%  of  its  length.  The  flutter  determinant  for  the  blade  is  repeatedly  solved  for 
multiple  elements  along  the  blade  model,  as  well  as  a  range  of  forward  velocities,  in  the 
conventional  manner  to  obtain  the  blade  flutter  speed  [Ref.  8].  UH-60  blade  parameters 
will  be  used  in  the  sample  problem  in  deference  to  the  extensive  database  now  available 
for  this  blade.  However,  for  the  demonstration  analysis  the  UH-60  blade  will  be 
modified  to  make  it  “flutter  susceptible”  by  moving  the  chord-wise  position  of  the  blade 
CG  aft  while  keeping  its  elastic  axis  at  the  quarter  chord.  This  also  introduces  flap 
torsion  coupling,  a  desirable  feature  for  a  sample  problem  of  this  type. 

For  the  numerical  example,  the  blade  will  be  divided  into  radial  segments  in  the 
usual  manner  with  unsteady  aerodynamics  applied  to  the  blade  at  each  panel  point.  The 
number  of  elements  chosen  has  some  effect  on  the  accuracy  of  the  solution,  since  the 
variation  of  velocity  along  the  span  must  be  adequately  incorporated.  Through  trial  and 
error,  it  has  been  found  that  the  change  in  flutter  point  as  the  number  of  elements  exceeds 
50  is  insignificant  (<1%),  while  using  only  5  or  10  elements  can  give  variations  as  great 
as  10%. 

Physically,  the  method  of  solution  in  this  thesis  is  equivalent  to  locking  the  blade 
at  the  90-degree  azimuth  position  and  solving  the  flutter  problem,  similar  to  a  fixed  wing 
case  with  Theodorsen  lift  deficiency  values,  yet  allowing  radial  velocity,  and  thus, 
reduced  frequency,  to  vary  with  span  as  in  the  case  of  the  tangential  velocity  of  a  rotor 

o 

blade  in  forward  flight,  \\f  =  90  . 
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The  average  characteristics  of  the  H-60  rotor  blade  are  given  in  Table  1  [Ref.  15]. 
Scanlan  and  Rosenbaum  [Ref.  8]  state  that  any  mode  with  a  natural  frequency  greater 
than  1 .2  times  the  frequency  of  the  first  torsion  mode  will  not  have  a  significant  effect  on 
the  flutter  point.  The  third  bending  mode  (>7.5f2)  and  the  second  torsion  mode  (>1  IQ) 
are  both  significantly  outside  this  range. 


Ta 


ble  1.  UH-60  Blade  Characteristics 


PARAMETER 

VALUE 

UNITS 

M 

0.00164 

Ib-s7irf 

C 

20.76 

in 

■« 

0.037 

Ib-sAn/in 

R 

322 

in 

E 

15 

in 

EA 

25%  chord 

CG 

Variable 

Q. 

27.02 

rad/s 

Q.R 

725 

ft/s 

Hence,  only  the  first  two  flap-wise  bending  modes  with  natural  frequencies  of 


2.8f2  and  4.712  are  incorporated  in  the  flutter  analysis  together  with  one  torsion  mode, 


occurring  at  4.3f2.  The  rotating  natural  frequencies  of  the  uniform  beam  model  are 
compared  to  that  of  the  real  blade  in  Table  2.  Note  that  the  frequencies  in  each  case  are 
within  10%  of  the  real  values,  giving  confidence  that  the  uniform  beam  model  is 
sufficient  for  a  first  order  approximation. 

Table  2.  Comparison  of  Rotating  Natural  Frequencies 


MODE 

BLADE  MODEL 

UH-60  BLADE 

1st  Bending 

72.75 

75.93 

2nd  Bending 

130.47 

140.23 

1st  Torsion 

128.84 

116.1 

With  the  mode  shapes  and  natural  frequencies  calculated,  the  flutter  analysis 
portion  of  the  program  yielded  results  for  a  range  of  velocities  for  each  CG  location 
specified,  from  the  quarter  chord  aft  to  the  %  chord.  After  this  analysis  was  conducted,  it 
was  done  again,  this  time  holding  the  forward  velocity  at  zero  and  using  the  same  tip 
velocity  range  as  the  forward  flight  case.  This  effort  simulates  whirl-stand  over-speed 
tests.  Figures  8  and  9  show  the  results  of  ro  vs.  Ufl  and  g  vs.  Ufl  for  both  the  forward 
velocity  and  over-speed  cases  with  zero  CG  offset. 
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Figure  8.  CG  Offset  Zero,  Forward  Flight 
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Figure  9.  CG  Offset  Zero,  Over-Speed 


28 


Note  that  the  blade  model  is  stable  in  both  cases,  with  the  natural  frequencies 
remaining  relatively  constant,  and  the  damping  coefficients  negative  throughout  the 
velocity  range.  Incidentally,  the  velocity  range  plotted  is  equivalent  to  the  forward  flight 
airspeeds  of  0  -  220kts  (approximately),  which  is  the  published  flight  envelope  for  the 
UH-60,  plus  40kts. 

Figures  10  through  21  show  the  progression  of  results,  in  the  same  fonnat,  as  the 
CG  is  moved  to  the  mid-chord  position  and  subsequently  aft  to  the  3A  chord  position. 
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Figure  10.  CG  Offset  0.5b,  Forward  Flight 
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Figure  11.  CG  Offset  0.5b,  Over-Speed 
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VELOCITY  (ft/s) 


2nd  Bending  1st  Torsion  —  1st  Bending 


VELOCITY  (ft/s) 


2nd  Bending  1st  Torsion  —  1st  Bending 
Figure  12.  CG  Offset  0.75b,  Forward  Flight 
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VELOCITY  (ft/s) 


2nd  Bending  1st  Torsion  —1st  Bending 


VELOCITY  (ft/s) 


2nd  Bending  1st  Torsion  —1st  Bending 
Figure  13.  CG  Offset  0.75b,  Over-Speed 
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Figure  14.  CG  Offset  0.85b,  Forward  Flight 
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VELOCITY  (ft/s) 


2nd  Bending  1st  Torsion  —1st  Bending 


Figure  15.  CG  Offset  0.85b,  Over-Speed 
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Figure  16.  CG  Offset  0.9b,  Forward  Flight 
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VELOCITY  (ft/s) 


2nd  Bending  1st  Torsion  —1st  Bending 


Figure  17.  CG  Offset  0.9b,  Over-Speed 
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VELOCITY  (ft/s) 


2nd  Bending  1st  Torsion  —  1st  Bending 


VELOCITY  (ft/s) 


2nd  Bending  1st  Torsion  —  1st  Bending 
Figure  18.  CG  Offset  0.95b,  Forward  Flight 
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VELOCITY  (ft/s) 


2nd  Bending  1st  Torsion  —1st  Bending 


Figure  19.  CG  Offset  0.95b,  Over-Speed 
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VELOCITY  (ft/s) 


2nd  Bending  1st  Torsion  —  1st  Bending 


VELOCITY  (ft/s) 


2nd  Bending  1st  Torsion  —  1st  Bending 
Figure  20.  CG  Offset  1.0b,  Forward  Flight 
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VELOCITY  (ft/s) 


2nd  Bending  1st  Torsion 


1st  Bending 


VELOCITY  (ft/s) 


2nd  Bending  1st  Torsion 


1st  Bending 


Figure  21.  CG  Offset  1.0b,  Over-Speed 


41 


B.  DISCUSSION  OF  RESULTS  FOR  UH-60  BLADE  EXAMPLE 


From  Figures  10  and  11,  it  can  be  seen  that  the  first  two  coupled  inodes, 
corresponding  to  the  first  bending  and  torsion  modes,  begin  to  converge  in  frequency, 
causing  one  of  the  modes  to  become  more  stable  while  the  other  becomes  less  stable. 
The  movement  of  the  damping  coefficient  curves,  show  this  stability,  or  instability:  larger 
negative  values  equals  increased  stability,  and  vice-versa.  The  first  flutter  point  appears 
for  the  over-speed  test  at  1017 ft/ s  tip  speed,  with  the  CG  at  0.75b,  or  62.5%  chord.  The 
forward  flight  curve  is  stable  at  this  point.  The  forward  flight  condition  remains  stable 
until  the  CG  is  moved  to  the  0.9b  position,  or  70%  chord.  The  tip  speed  for  this  case  is 
1054ft/s,  corresponding  to  a  forward  flight  velocity  of  194kts.  Note  that  the  flutter  point 
for  each  test  occurs  at  a  lower  tip  speed,  as  the  CG  is  moved  further  aft.  Table  3  provides 
a  summary  of  the  flutter  points,  for  both  the  maximum  forward  flight  velocity  and 
maximum  over-speed.  The  forward  flight  column  lists  the  tip  speed  and  forward  velocity 
of  the  flutter  point,  and  the  over-speed  column  lists  the  flutter  point  tip  speed. 


Table  3.  Flutter  Points 


CG  Offset 

CG  (%  chord) 

Forward  Flight 

Over-speed 

0.75b 

62.5% 

- 

1017fps 

0.85b 

67.5% 

- 

960fps 

0.90b 

70.0% 

1054fps/  194kts 

937fps 

0.95b 

72.5% 

981fps  /  15  lkts 

917fps 

1.00b 

75.0% 

933fps  /  123kts 

893fps 
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IV.  SUMMARY 


A.  CONCLUSIONS 

This  thesis  has  presented  a  fully  coupled,  quasi-3D  analysis  of  rotor  blade  flutter 
that  can  accommodate  forward  flight  conditions.  The  rotor  blade  is  modeled  as  a  uniform 
beam,  taking  the  average  characteristics  of  a  real  blade  between  20%  and  90%  of  its 
length.  Applying  Rayleigh’s  method,  the  first  few  bending  and  torsion  normal  mode 
shapes  and  natural  frequencies  are  determined,  and  then  adjusted  for  the  rotating  case. 
With  this  data,  force  and  moment  equations  of  motion  are  developed  using  LaGrange’s 
equation  along  with  a  normal  mode  analysis.  Theodorsen’s  lift  deficiency  function  is 
calculated  over  a  range  of  forward  velocities,  or  reduced  frequencies,  for  a  specified 
number  of  elements  along  the  blade  model.  Incorporating  these  coefficients  into  the 
equations  of  motion,  a  square  matrix  is  generated  from  which  complex  eigenvalues  can 
be  derived.  These  eigenvalues  provide  the  aeroelastic  natural  frequencies  and  damping 
coefficients  for  each  coupled  mode.  The  forward  velocity  at  which  one  of  the  modes 
produces  a  damping  coefficient  of  zero  gives  the  value  of  reduced  frequency  for  the 
flutter  point.  The  resulting  forward  speed  and  blade  tip  speed  can  then  be  determined. 

The  results  of  this  analysis  appear  to  validate  the  methods  used.  The  UH-60 
model  blade  remains  stable  within  a  feasible  forward  velocity  range  (0  to  180kts)  until 
the  CG  is  moved  aft  of  70%  chord.  As  the  CG  is  moved  aft,  two  of  the  coupled  modes, 
corresponding  to  the  first  torsion  and  first  bending  modes,  converge  on  a  single 
frequency,  approximately  91.8rad/s.  As  this  convergence  approaches,  the  mode 
corresponding  to  the  first  torsion  mode  becomes  decreasingly  stable,  while  the  other 
becomes  increasingly  stable.  Once  the  modes  converge  sufficiently,  due  to  CG  location, 
the  coupled  first  torsion  mode  rapidly  becomes  unstable  in  the  characteristic  ‘S-curve’ 
manner.  As  a  final  indication  of  the  validity  of  this  analysis,  the  flutter  point  occurs  at 
decreasing  forward  speeds  as  the  CG  is  moved  even  further  aft  than  the  70%  chord  point. 
The  over-speed  method  of  approximating  forward  flight  flutter  conditions,  in  contrast  to 
the  forward  flight  analysis  presented  above,  seems  to  provide  more  conservative  results. 
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The  method  presented  in  this  thesis  appears  to  be  a  valid  and  useful  cost-  and 
weight-saving  tool  for  the  rotor  blade  designer.  This  relatively  simple  approach  can 
provide  vital  infonnation  for  designing  whirl-stand  as  well  as  flight  tests  of  newly 
developed  rotor  blades,  without  being  excessively  conservative.  The  weight  savings 
provided  by  allowing  CG  repositioning  can  be  used  to  optimally  weight  the  blade  for 
better  autorotative  performance,  or  decreased  design/weight/cost  expenditures  on  blade 
damping.  Altogether,  further  development  of  pre-existing  flutter  theory  and  unification 
with  the  power  of  modern  computers  allows  advanced  prediction  without  extensive  cost. 
B.  RECOMMENDATIONS 

Some  of  the  simplifications  made  are  not  necessary,  given  the  power  and  speed  of 
modern  computers.  It  should  be  noted  that  the  simplifications  of  uniform  blade  and 
Yntema’s  method  were  used  here  to  make  possible  rapid  parametric  studies.  For  detailed 
analysis  of  an  actual  rotor  blade,  the  blade  frequencies  and  mode  shapes  can  be  more 
accurately  determined  by  a  transfer  matrix  approach,  such  as  the  Mykelstad  and/or 
Holtzer  methods.  In  addition,  the  incorporation  of  more  advanced  lift  deficiency 
functions,  such  as  those  developed  by  Loewy,  and  Shipman  and  Wood,  which  take  into 
account  the  contributions  of  an  unsteady  wake,  could  further  increase  the  accuracy  of 
flutter  prediction. 

Another  area  to  consider  is  that  compressibility  effects  are  not  taken  into  account 
in  spite  of  the  fact  that,  even  in  a  hover,  the  tip  speed  of  most  rotor  blades  is  well  into  the 
compressible  region.  At  high  forward  velocities,  the  tip  speed  often  approaches  Mach  1 , 
and  is  likely  to  encounter  transonic  effects  as  well.  The  incorporation  of  compressibility 
effects  promises  to  be  challenging,  but  is  likely  to  increase  the  accuracy  of  prediction. 
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APPENDIX  A.  3D  BLADE  FLUTTER  PROGRAM 


A.  PROGRAM  WALK-THROUGH 

A  MATLAB™  program  was  written  to  conduct  the  flutter  analysis  using  a 
uniform  beam  model  of  a  real  blade.  MATLAB™  provides  the  user  with  the  ability  to 
make  calculations  simultaneously  for  large  multi-dimensional  matrices,  making  it 
appropriate  for  this  subject  matter.  The  steps  of  the  program,  along  with  an  explanation 
of  the  methodology,  follows: 

1.  Step  1  -  User  Input 

The  program  begins  by  clearing  the  workspace  variables.  Then  the  user  must 
choose  the  number  of  modes  (3  or  5),  the  number  of  blade  elements  and  the  altitude  (for 
air  density  variation).  The  separation  between  the  CG  and  EA  (positive  for  CG  aft  of 
EA)  and  the  estimated  forward  velocity  range  must  also  be  chosen.  For  the  sample 
problem,  a  three-mode  analysis  was  conducted  over  100  blade  elements  at  sea  level. 

2.  Step  2  -  Establish  Blade  Model  Properties 

The  UH-60  blade  was  modeled  as  a  unifonn  beam.  The  values  chosen  to 
represent  the  blade  are  the  average  values  between  20%  and  90%  of  the  blade’s  length. 
Properties  needed  for  this  analysis  include  the  mass  per  unit  length  (p),  Young’s  Modulus 
(E),  the  flap-wise  (EIXX)  and  polar  (Ia)  mass  moments  of  inertia,  and  torsional  stiffness 
(GJ).  Throughout  the  program,  the  units  of  important  quantities  are  denoted  for 
clarification.  In  addition,  this  step  establishes  important  aspects  of  the  UH-60  such  as 
rotor  speed  and  blade  offset,  as  well  as  the  location  of  the  elastic  axis  (in  percent  semi¬ 
chord  from  the  mid-chord). 

3.  Step  3  -  Calculate  Mode  Shapes  and  Natural  Frequencies 

This  step  establishes  the  elastic  and  inertial  properties  of  the  uniform  beam  and 
displays  them,  for  comparison  to  real  blade  properties.  First,  the  constants  for  the 
analysis  are  set.  The  references  for  these  constants  are  noted  in  the  program  itself,  and 
explained  in  the  section  on  theoretical  development.  The  beam  is  then  divided  into  the 
user-specified  number  of  elements,  and  the  radii  to  the  midpoints  of  the  elements  are 
calculated.  The  midpoints  are  the  locations  where  the  2-D  calculations  will  be  made. 
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The  mode  shapes  for  bending  and  torsion  are  then  calculated  and  normalized  to  a 
maximum  deflection  of  one  and  a  positive  initial  slope.  These  mode  shapes  are  shown  in 
Fig.  22.  A  small  ‘f  denotes  a  bending  mode  and  a  capital  ‘F’  denotes  a  torsion  mode. 


Figure  22.  Normal  Bending  and  Torsion  Mode  Shapes 

The  generalized  mass  for  each  mode  is  now  calculated,  using  the  mode  shapes 
just  obtained.  The  integral  for  the  generalized  masses  is  approximated  by  a  summation  as 
in  the  following  example: 

M2  =  J0m(x)[f2  (x)]” dx  =  &df2  (n) 

n=l 

where  p  is  the  number  of  elements  and  d  is  the  length  of  each  element.  Finally,  the  non¬ 
rotating  and  rotating  natural  frequencies  are  calculated  with  the  equations  above  and 
displayed  for  comparison  with  values  from  the  real  blade. 
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4.  Step  4  -  Construct  Matrices  for  Velocity  Range 

The  fourth  step  is  the  heart  of  the  program.  This  section  begins  the  loop  that 
calculates  the  matrix  and  its  eigenvalues  for  each  forward  velocity.  The  range  of 
velocities  is  set  by  the  user’s  estimate  and  divided  into  300  equally  spaced  velocities. 
This  number  seems  to  provide  good  resolution  throughout  the  range  of  calculations. 
Next,  a  set  of  reduced  frequencies  is  calculated  using  the  semi-chord  and  the  first  torsion 
mode  rotating  natural  frequency: 

k  _  bo^ 

"P  V„  +  £>R 

This  frequency  is  used  as  an  approximation  to  the  actual  coupled  frequency  that  is 
expected  to  become  unstable  and  is  only  used  as  a  means  to  obtain  a  reasonable  range  of 
reduced  frequencies. 

Following  this,  the  static  moment  coupling  terms  must  be  calculated.  Again,  the 
integral  is  approximated  with  a  summation  as  follows: 

Sa12  =  JA  (X)fl  (X)F2  (X)dx  =  XSa  (n)fl  (n)F2  (n)d 

n=l 

where  p  is  the  user-specified  number  of  elements,  and  d  is  the  length  of  each  element. 

Now  the  program  is  set  to  begin  iterative  calculations.  For  each  step  of  forward 
velocity,  Theodorsen  coefficients  are  calculated  for  every  element  along  the  blade, 
according  to  the  differing  local  velocities,  and  thus,  reduced  frequencies.  Arrays  are 
created  of  reduced  frequencies,  lift  deficiency  functions,  and  coefficients,  corresponding 
to  the  array  of  element  midpoints.  Then,  the  matrix  of  integrals  in  equation  set  (5)  is 
constructed  as  a  set  of  summations  in  the  same  manner  as  previously  noted.  The  program 
variable  for  this  matrix  is  ‘A’.  The  ‘A’  matrix  is  then  modified  to  incorporate  the 
dimensional  quantities  of  air  density,  generalized  mass,  and  coupled  static  moment,  as 
well  as  the  factor  7t.  This  matrix,  denoted  ‘Adim,’  is  analogous  to  the  matrix  given  in 
equation  (6).  Finally,  the  ‘Adim’  matrix  is  ‘re-dimensionalized’  with  the  generalized 
masses  and  squared  natural  frequency  ratios  in  order  to  generate  eigenvalues  of  the 
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desired  form.  The  resulting  matrix,  ‘Abar’  in  the  program,  is  the  final  five  by  five, 
generalized  forces  and  moments  given  in  equation  set  (7). 

5.  Step  5  -  Calculate  Eigenvalues 

The  user-specified  number  of  modes  must  be  accounted  for  before  the 
eigenvalues  can  be  calculated.  If  the  user  picks  five  modes,  the  full  matrix  will  be  used, 
but  this  currently  yields  a  result  that  has  one  eigenvalues  with  a  negative  real  part,  which 
corresponds  to  an  imaginary  frequency  that  is  not  physically  permissible.  Therefore,  the 
three-mode  matrix  must  be  used.  Luckily,  no  further  calculations  need  be  made  to  enact 
this  solution.  A  new  matrix  is  fonned  from  the  5  by  5  matrix  by  simply  eliminating  all 
tenns  with  the  subscripts  3  and  5,  which  correspond  to  the  third  bending  mode  and  the 
second  torsion  mode,  creating  the  appropriate  3  by  3  coupled  matrix.  Conveniently 
enough,  it  is  noted  in  Scanlan  and  Rosenbaum  [Ref.  7]  that  any  mode  with  a  natural 
frequency  greater  than  1 .2  times  the  frequency  of  the  first  torsion  mode  will  not  have  a 
significant  effect  on  the  flutter  point.  Both  of  the  eliminated  modes  fall  into  this 
category. 

Finally,  the  eigenvalues  are  calculated  using  MATLAB™’s  ‘eig’  command.  This 
command  is  executed  without  balancing  in  order  to  eliminate  the  possibility  that  a  small 
value  off  the  diagonal  will  be  inordinately  skewed  in  the  balancing  process.  The 
eigenvalues  are  then  sorted  according  to  their  real  part  and  broken  down  into  coupled 
frequencies  and  damping  coefficients  as  discussed  above. 

6.  Step  6  -  Display  Results  and  Save  Data 

The  last  step  is  to  save  and  display  the  data  in  a  useful  manner.  Recall  that  the 
values  of  reduced  frequency  were  obtained  using  the  natural  frequency  of  the  first  torsion 
mode.  However,  the  coupled  frequency  corresponding  to  the  first  torsion  mode  varies 
with  velocity  in  a  non-linear  fashion.  To  obtain  the  actual  velocity  for  each  point 
calculated,  the  range  of  reduced  frequencies  must  be  divided  into  the  semi-chord  and  the 
range  of  coupled  (1st  torsion  mode)  frequencies. 


Note  that  the  division  is  conducted  term  by  term,  as  an  array. 
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The  frequencies  and  damping  ratios  are  then  plotted  versus  this  velocity 
such  that  the  user  can  view  the  flutter  point  (damping  coefficient  ‘g’  =  0),  if  any.  The 
data  is  also  saved  to  a  file  for  incorporation  into  a  spreadsheet  for  further  data  analysis,  as 
desired. 


B.  PROGRAM  LISTING 


%  Program  to  find  natural  static  &  rotating 
%  mode  shapes  &  frequencies  of  a  unifonn  beam 
%  Used  to  model  the  UH-60  rotor  blade 
%  (no  aerodynamic  effects  taken  into  account) 

%  Following  frequency  determination,  all  effects 
%  are  taken  into  account  to  determine 
%  the  flutter  point  (g=0  for  one  mode) 
clear 

mode=3;  %  input('5x5  (5)  or  3x3  (3)?  ’); 

p=100;  %  input(’How  many  elements?  ’); 

alt=0;  %  input(’What  altitude?  (0,  5k,  10k)  ’); 

cgea=input(’CG  -  EA  separation?  ’); 

mn=l ;  %  input(’Vfwd  minimum?  ’); 

mx=90 1 ;  %  input(’V fwd  maximum?  ’); 


R=322;  offset=15;  chord=20.76; 
L=R-offset;  e=offset/L; 

0=27.02; 

mu=0.00164167; 

E=10000000; 

EIxx=2.27859e7;  EIyy=7.71261e8; 

I=(EIxx+EIyy)/E; 

la=0. 037006386; 

GJ=2.467909e7; 


%  inches 
%  inches 

%  Rotor  Speed  rad/s 
%  lbsA2/inA2 
%  lb/inA2 
%  lbinA2 
%  inA4 

%  lbinsA2/in  polar  mass  MOI  about  e.a. 
%  lbinA2 


a=-0.5; 
b=chord/2; 
if  alt==l  0 

rho=0 . 00 1 75  5 6/(  1 2 A4) ; 
elseif  alt=5 
rho=0.0020482/(  12A4); 
else 

rho=0.0023769/(  12A4); 
end 


%  non-dim  elastic  axis  location  measured  from  midchord 
%  Half  chord  length  in  inches 

%  10000  ft  %  lbsA2/inA4  density  of  air 

%  5000  ft 

%  Sea  Level 


%  for  rho  =  0  AND  eg  -  ea  =  0  :  the  natural  frequencies  are  returned 
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asubn=[  15.4 1820562  49.96486209  104.2476966];  %  Hartog 

K0=[6.38  17.63  35.05];  Kl=[9.18  26.02  52.2];  %  Yntema 

BnR=[3. 926602  7.068583  10.21018];  %  Young/Felgar  Free-Supported 

An=[l. 000777  1.000001  1];  %  Young/Felgar  Free-Supported 

d=R/p;  %  length  in  inches  of  one  element 

Salf=mu*d*cgea*b;  %  Calculate  the  static  moment  of  each  element 

r(l)=0;  %  r  is  the  vector  of  element  edges  in  percent  R 

for  q=l:p 
r(q+l)=q/p; 
end 

for  n=T:p 

mid(n)=(r(n+ 1  )+r(n))/2 ;  %  midpoints  of  elements 

end 

u=mid'*BnR;  %  This  section  calculates  the  bending  mode  shapes 

for  n=l  ;3  %  According  to  Young  and  Felgar  Supported-Free 

for  m=l:p 
h=u(m,n); 

f(p+l-m,n)=(cosh(h)+cos(h)-An(n)*(sinh(h)-l-sin(h)))/2; 

%  divide  by  two  to  normalize  max  deflection  to  one 
end 
end 

for  n=l:3 
if  f(  1  ,n)<0 
h=-l*f(:,n); 
f(:,n)=h; 
end 
end 

for  n=l  :3  %  This  section  calculates  the  torsional  mode  shapes 

for  m=l:p  %  According  to  Hartog  (already  normalized  to  1) 

F(m,n)=sin((n-.5)*pi*mid(m)); 
end 
end 

figure(  1)  %  Plot  of  mode  shapes 

plot(mid,f(:,l),mid,f(:,2),mid,f(:,3),niid,F(:,l),mid,F(:,2)) 
legend('f(  l)','f(2)','f(3)','F(  1)','F(2)') 

f$qr=f.*f;  %  Calculate  Generalized  Mass  for  bending  and  torsion  modes 
FSqr=F.*F; 
fMsqr=mu*  d*  IBqr; 

FMsqr=d*Ia*FSqr; 


%  This  section  normalizes  mode  shape 
%  to  all  begin  with  a  positive  slope 
%  for  display  purposes 
%  (does  not  affect  calculations) 
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Gmass=sum(fMsqr); 

GTmass=sum(FMsqr); 

oN=asubn*sqrt(EIxx/(mu*RA4));  %  Calculate  natural  non-rotating  freqs 

oTN=[0.5  1.5  2.5]*pi*sqrt(GJ/(Ia*RA2)); 

oR=sqrt(oN.*oN+(K0+Kl*e).*OA2);  %  Calculate  natural  rotating  freqs 

oTR=(0.A2+oTN.A2).A(  1/2); 

W=[Gmass;  oN;  oR]’ 

X=[GTmass;  oTN;  oTR]’ 

Vfwd=linspace(mn,mx,300)*  12; 

ktip=b  *  oTR(  1 ) ./( Vfwd+O  *  R) ; 

for  q=l:length(ktip) 

H 1  (q)=besselh(  1 ,2,ktip(q)); 

H0(q)=besselh(0,2,ktip(q)); 

C(q)=H  1  (q)/ (H 1  (q)+i*H0(q)); 

RE(q)=real(C(q)); 

IM(q)=imag(C(q)); 
end 

figure(2)  %  Plot  Theodorsen  lift  deficiency  function  for  our  range 
semilogx(ktip,RE,ktip,-IM) 
legend('FV-G’); 

for  n=l  :3  %  Calculate  static  moments  for  mode  coupling 

for  m=l:2 

Sa(n,m)=sum(SalPf(:,n).*F(:,m)); 
end 
end 


%  Display  bending  results 
%  Display  torsion  results 

%  Set  range  of  airspeeds  for  helicopter 
%  Define  reduced  frequency  for  tip 

%  validate  Theodorsen  calculations 
%  by  comparing  w/  F  and  -G  plots 
%  from  prior  works 


pirho=pi*rho; 

oal=oTR(l); 

for  q=l  :length(Vfwd)  %  Run  entire  blade  for  range  of  values  k  at  blade  tip 

for  t=l  :length(mid)  %  Calc  Theodorsen  coeff  s  for  each  element  of  blade 

K(t)=(b  *  oa  1 ) ./( V  fwd(q)+(0  *R*mid(p+ 1  -t))); 

HI  (t)=besselh(  1 ,2,K(t)); 

H0(t)=besselh(0,2,K(t)); 

C  (t)=H  1  (t)/(H  1  (t)+i  *H0(t)) ; 

Lh(t)=(  1  -2  *  i  *  C(t)/K(t)) ; 

La(t)=(0.5-2*i*(0.5+(l-i/K(t))*C(t))/K(t)); 

Mh(t)=(0.5); 

Ma(t)=(3/8-i/K(t)); 
end 


%  Theodorsen’ s  function 
%  Lift  -  plunge 
%  Lift  -  pitch 
%  Moment  -  plunge 
%  Moment  -  pitch 
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%  %  %  %  %  Calculate  A  matrix  (Generalized  Force  Matrix) 
A(l,l)=sum(d*Lh(:)*bA2.*f(:,l).A2);  %  inA3 

A(  1 ,2)=sum(d* Lh( : )  *  b  A2 .  *  f( : ,  1 ) .  *  f( :  ,2));  %  inA3 

A(l,3)=sum(d*Lh(:)*bA2.*f(:,l).*f(:,3));  %  inA3 

A(l,4)=sum(d*(La(:)-(0.5+a)*Lh(:))*bA3.*f(:,l).*F(:,l)); 
A(l,5)=sum(d*(La(:)-(0.5+a)*Lh(:))*bA3.*f(:,l).*F(:,2)); 
A(2,1)=A(1,2);  %  inA3 

A(2,2)=sum(d*Lh(:)*bA2.*f(:,2).A2);  %  inA3 

A(2,3)=sum(d*Lh(:)*bA2.*f(:,3).*f(:,2));  %  inA3 

A(2,4)=sum(d*(La(:)-(0.5+a)*Lh(:))*bA3.*f(:,2).*F(:,l)); 
A(2,5)=sum(d*(La(:)-(0.5+a)*Lh(:))*bA3.*f(:,2).*F(:,2)); 
A(3,1)=A(1,3);  A(3,2)=A(2,3);  %  inA3 

A(3,3)=sum(d*Lh(:)*bA2.*f(:,3).A2);  %  inA3 

A(3,4)=sum(d*(La(:)-(0.5+a)*Lh(:))*bA3.*f(:,3).*F(:,l)); 
A(3,5)=sum(d*(La(:)-(0.5+a)*Lh(:))*bA3.*f(:,3).*F(:,2)); 
A(4,l)=sum(d*(Mh(:)-(0.5+a)*Lh(:))*bA3.*f(:,l).*F(:,l)); 
A(4,2)=sum(d*(Mh(:)-(0.5+a)*Lh(:))*bA3.*f(:,2).*F(:,l)); 
A(4,3)=sum(d*(Mh(:)-(0.5+a)*Lh(:))*bA3.*f(:,3).*F(:,l)); 
A(4,4)=sum(d*(Ma(:)-(0.5+a)*(La(:)+Mh(:))+(0.5+a)A2*Lh(:)) 
*bA4.*F(:,l).A2); 

A(4,5)=sum(d*(Ma(:)-(0.5+a)*(La(:)+Mh(:))+(0.5+a)A2*Lh(:)) 

*bA4.*F(:,l).*F(:,2)); 

A(5,l)=sum(d*(Mh(:)-(0.5+a)*Lh(:))*bA3.*f(:,l).*F(:,2)); 

A(5,2)=sum(d*(Mh(:)-(0.5+a)*Lh(:))*bA3.*f(:,2).*F(:,2)); 

A(5,3)=sum(d*(Mh(:)-(0.5+a)*Lh(:))*bA3.*f(:,2).*F(:,2)); 

A(5,4)=A(4,5); 

A(5,5)=sum(d*(Ma(:)-(0.5+a)*(La(:)+Mh(:))+(0.5+a)A2*Lh(:)) 

*bA4.*F(:,2).A2); 


%  inA4 
%  inA4 


%  inA4 
%  inA4 


%  inA4 
%  inA4 
%  inA4 
%  inA4 
%  inA4 

%  inA5 

%  inA5 
%  inA4 
%  inA4 
%  inA4 
%  inA5 

%  inA5 


%  %  %  %  %  Incorporate  pi,  rho,  generalized  masses  and  MOI  correction 


Adim(  1 , 1  )=pirho  *  A(  1 , 1  )+Gmass(  1 ) ; 

% 

lbsA2/in 

Adim(  1 ,2)=pirho*  A(  1 ,2); 

% 

lbsA2/in 

Adim(  1 ,3)=pirho*  A(  1 ,3); 

% 

lbsA2/in 

Adim(  1 ,4)=pirho*  A(  1 ,4)+Sa(  1,1); 

% 

lbsA2 

Adim(  1 ,5)=pirho*  A(  1 ,5)+Sa(  1 ,2); 

% 

lbsA2 

Adim(2,l)=pirho*A(2,l); 

% 

lbsA2/in 

Adim(2,2)=pirho*A(2,2)+Gmass(2); 

% 

lbsA2/in 

Adim(2,3)=pirho*A(2,3); 

% 

lbsA2/in 

Adim(2,4)=pirho*A(2,4)+Sa(2,l); 

% 

lbsA2 

Adim(2,5)=pirho*A(2,5)+Sa(2,2); 

% 

lbsA2 

Adim(3 , 1  )=pirho  *  A(3 , 1 ) ; 

% 

lbsA2/in 

Adim(3 ,2)=pirho  *  A(3 ,2); 

% 

lbsA2/in 

Adim(3 ,3)=pirho*  A(3 ,3  )+Gmass(3); 

% 

lbsA2/in 

Adim(3 ,4)=pirho*  A(3 ,4)+Sa(3 , 1 ); 

% 

lbsA2 

Adim(3,5)=pirho*A(3,5)+Sa(3,2); 

% 

lbsA2 

Adim(4, 1  )=pirho*  A(4, 1  )+Sa(  1,1); 

% 

lbsA2 
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Adim(4,2)=pirho*A(4,2)+Sa(2,l); 

% 

lbsA2 

Adim(4,3)=pirho*A(4,3)+Sa(3,l); 

% 

lbsA2 

Adim(4,4)=pirho  *  A(4,4)+GT  mass(  1 ) ; 

% 

lbinsA2 

Adim(4,5)=pirho*A(4,5); 

% 

lbinsA2 

Adim(5 , 1  )=pirho*  A(5 , 1  )+Sa(  1 ,2); 

% 

lbsA2 

Adim(5 ,2)=pirho*  A(5 ,2)+Sa(2,2); 

% 

lbsA2 

Adim(5,3)=pirho*A(5,3)+Sa(3,2); 

% 

lbsA2 

Adim(5 ,4)=pirho  *  A(5 ,4); 

% 

lbinsA2 

Adim(5,5)=pirho*A(5,5)+GTmass(2); 

% 

lbinsA2 

%  %  %  %  %  Manipulate  A  matrix  to  be  eigenvalue  friendly 
Abar(  1 , 1  )=(Adim(  1 , 1  )/Gmass(  1 ))  *  (oTR(  1  )/oR(  1  ))A2; 

Abar(  1 ,2)=(Adim(  1 ,2)/Gmass(l))*(oTR(  l)/oR(  1))A2; 

Abar(  1 ,3)=(Adim(  1 ,3)/Gmass(l))*(oTR(  l)/oR(  1))A2; 

Abar(  1 ,4)=(Adim(  1 ,4)/ Gmass(  1 ))  *  (oTR(  1 )/ oR(  1  ))A2;  %  in 

Abar(l,5)=(Adim(l,5)/Gmass(l))*(oTR(l)/oR(l))A2;  %  in 

Abar(2 , 1  )=(Adim(2, 1 )/ Gmass(2))  *  (oTR(  1 )/ oR(2))A2; 
Abar(2,2)=(Adim(2,2)/Gmass(2))*(oTR(l)/oR(2))A2; 
Abar(2,3)=(Adim(2,3)/Gmass(2))*(oTR(l)/oR(2))A2; 
Abar(2,4)=(Adim(2,4)/Gmass(2))*(oTR(l)/oR(2))A2;  %  in 

Abar(2,5)=(Adim(2,5)/Gmass(2))*(oTR(l)/oR(2))A2;  %  in 

Abar(3 , 1  )=(Adim(3 , 1  )/Gmass(3))  *  (oTR(  1  )/oR(3))A2; 
Abar(3,2)=(Adim(3,2)/Gmass(3))*(oTR(l)/oR(3))A2; 
Abar(3,3)“(Adim(3,3)/Gmass(3))Hs(oTR(l)/oR(3))A2; 
Abar(3,4)=(Adim(3,4)/Gmass(3))*(oTR(l)/oR(3))A2;  %  in 

Abar(3,5)=(Adim(3,5)/Gmass(3))*(oTR(l)/oR(3))A2;  %  in 

Abar(4,l)=Adim(4,l)/GTmass(l);  %  /in 

Abar(4,2)=Adim(4,2)/GTmass(l);  %  /in 

Abar(4,3)=Adim(4,3)/GTmass(l);  %  /in 

Abar(4,4)=Adim(4,4)/GTmass(l); 
Abar(4,5)=Adim(4,5)/GTmass(l); 

Abar(5 , 1  )=(Adim(5 , 1 )/ GTmass(2))*  (oTR(  1  )/oTR(2))A2 ;  %  /in 
Abar(5,2)=(Adim(5,2)/GTmass(2))*(oTR(l)/oTR(2))A2;  %  /in 
Abar(5,3)=(Adim(5,3)/GTmass(2))*(oTR(l)/oTR(2))A2;  %  /in 
Abar(5,4)=(Adim(5,4)/GTmass(2))*(oTR(l)/oTR(2))A2; 
Abar(5,5)=(Adim(5,5)/GTmass(2))Hs(oTR(l)/oTR(2))A2; 

%  Define  Matrix  for  3  or  5  mode  analysis 
if  mode==3 

T(l,l)=Abar(l,l);  T(l,2)=Abar(l,2);  T(l,3)=Abar(l,4); 
T(2,l)=Abar(2,l);  T(2,2)=Abar(2,2);  T(2,3)=Abar(2,4); 
T(3,l)=Abar(4,l);  T(3,2)=Abar(4,2);  T(3,3)=Abar(4,4); 
else 

T=Abar; 

end 
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Z(:,q)=eig(T, ’nobalance');  %  Calculate  eigenvalues  w/o  conditioning 

[Zz(:,q),index]=sort(real(Z(:,q)));  %  sort  by  real  part  (frequency  part) 
Zl(:,q)=Z(index,q);  %  must  keep  real  &  imag  parts  together 

ZR( : ,  q)=real(Z  1 ( :  ,q)) ; 

ZI(:,q)=imag(Zl(:,q)); 

o( :  ,q)=real(oTR(  1 )  ./sqrt(ZR( :  ,q))); 

%  the  assumption  of  real  frequencies  is  implicit  p.340  S&R 
g(:,q)=ZI(:,q)./ZR(:,q);  %  damping  matrix 

end 

%  Generate  x-axis  values  for  plotting 
%  Only  need  to  calculate  values  for  the  unstable  mode 
%  Since  only  the  flutter  point  is  a  true  representation  of  reality 
x=l./ktip; 
if  mode==3 

V2=(o(2,q)*b/12).*x;  %  velocity  corresponding  to  actual  oTR(l)  in  fwd  flight 
else 

V2=(o(4,q)*b/12).*x;  %  velocity  corresponding  to  actual  oTR(l)  in  fwd  flight 
end 

%  Save  output  file  for  incorporation  into  spreadsheet  program  (better  plotting) 
spread=[g;  o;  V2]'; 

dlmwrite(’APlot  1 00.  txf, spread, '\f,  1,1); 

%  Display  graphs  of  mode  fireqs  &  damping  coeff  s  versus  tip  velocity  and  1/Ktip 

figure(3); 

subplot(2,l,l) 

plot(V2,g) 

axis([700, 1100,-0.5,0.1]) 
grid  on 

xlabel(’Tip  Velocity  (ft/s)');  ylabel(’damping  coefficient'); 

subplot(2,l,2) 

plot(V2,o) 

axis([700, 1100,50, 150]) 
grid  on 

xlabel(’Tip  Velocity  (ft/s)');  ylabel('frequency  of  oscillation'); 

figure(4) 

subplot(2,l,l) 

semilogx(x,g) 

grid  on 

xlabel(’l  /  Ktip');  ylabel(’damping  coefficient'); 

axis([min(x),max(x), -0.5,0. 1]) 

subplot(2,l,2) 

semilogx(x,o) 

grid  on 

xlabel(’l  /  Ktip');  ylabel('frequency  of  oscillation’); 
axis([min(x),max(x),50,150]) 
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